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Abstract. We minimize a discrete version of the fourth-order curvature based Landau free energy by 
extending Brakke's Surface Evolver. This model predicts spherical as well as non-spherical shapes with 
dimples, bumps and ridges to be the energy minimizers. Our results suggest that the buckling and faceting 
transitions, usually associated with crystalline matter, can also be an intrinsic property of non-crystalline 
membranes. 



1 Introduction 

The study of organized structures, like self-assembled mem- 
branes and vesicles, is an important part of soft matter 
physics that is relevant for chemistry and biology and for 
practical applications [112] . These systems can usually be 
described as two-dimensional surfaces, since their thick- 
ness is much smaller than the typical size. Thus, at the 
phenomenological level, the search for equilibrium shapes 
of such micro- and nano-structures is a problem of differ- 
ential geometry of surfaces. This view is certainly appro- 
priate for liquid membranes where the free energy is con- 
sidered as a functional of the local curvatures and does 
not depend on the in plane deformation tensor like for 
crystalline membranes [3j. 

A simple theory for elastic shells was proposed by So- 
phie Germain around 1810 [415] in which the energy is 
given by the form 



E(S) 



fI dS< " 



AH 2 + BK}, 



(1) 



where H is the mean curvature, K is the Gaussian cur- 
vature and the integral is over the surface S. In the case 
of soap films, the first term, proportional to the surface 
tension a, is the dominant contribution to the free energy. 
For fluid membranes, instead, a ~ because molecules 
can easily flow and adjust the total area of the mem- 
brane to the one corresponding to the best packing [3J. 
The second term in the integrand in Eq. ([T]) describes 
out-of-plane bending of elastic shells. According to the 
Gauss-Bonnet theorem [8] for compact surfaces, the last 
term JJ dS K = 2irx(S) is a topological invariant, with 
x(S) called the Euler— Poincare characteristic of the sur- 
face. For all surfaces topologically equivalent to spheres 
X = 2 and JJ dSK = 4tt. 

The squared mean curvature integral W — JJ dS H 2 is 
known in mathematics as Willmore functional [7J, in the 



theory of membranes as Helfrich free energy [5], and in 
string theory as Polyakov action [3J . The sphere turns out 
to give the absolute minimum of this functional, W = &tt, 
among all compact surfaces [TJ. This naturally explains 
why spherical shapes are common in the world of fluid 
membranes. More complicated equilibrium shapes than 
the spherical one, e.g. prolates, discocytes and stomato- 
cytes [9110] . can be obtained within the Willmore func- 
tional by adding the condition of constant volume. Com- 
plicated equilibrium shapes are also observed in systems, 
like red blood cells or viral capsids, where other degrees 
of freedom rather than bending become important [3J so 
that the functional W and thus Eq. ([I]) is not sufficient. 
These situations are usually described in terms of crys- 
talline membranes where these additional degrees of free- 
dom are described by the in-plane deformation tensor. 

In Ref . [llj , for symmetric bolaamphiphilic fluid mem- 
branes, it was found that the interactions of the hydrophilic 
tails of the bolaamphiphiles with molecules of the solvent 
as well as entropic terms may lead to a negative coefficient 
A in Eq. ([1]). Then, symmetry allowed higher order terms 
should be added to stabilize the free energy that other- 
wise would not be bounded from below leading to the free 
energy functional 



Ta = 



dS {-A^ + BK+d^+c^K+csK 2 }, (2) 



where the minus sign is written explicitly so that from 
now on A is positive. The interplay between these higher 
order terms and the Willmore functional with negative 
sign leads to spontaneous bending. This model was suc- 
cessfully applied to describe the experimental data on de- 
formations of spherical bolaamphiphilic vesicles in high 
magnetic fields [12] . Here we study the equilibrium shapes 
of liquid membranes with spontaneous bending based on 
minimization of the functional Eq. ([2]). 
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Without referring to any specific system, we minimize 
numerically J-4 under the constraint of constant surface 
(S = AttR 2 ) and preserving the topology of the sphere 
(X = 2). To this purpose we have supplemented the open- 
source software "Surface Evolver" [T3] with new subrou- 
tines for the calculation of fourth order terms. Increasing 
A, for a given set of Cj, we find a transition from a spher- 
ical shape to more complicated, dimpled, shapes when 
AR 2 /c\ ~ 2. This continuous transition is followed by 
a discontinuous one towards shapes of icosahedral sym- 
metry with ridges and facets, when AR 2 /c\ > 8. The 
presence of the coupling term between the mean curva- 
ture and the Gaussian curvature (c 2 7^ 0) results in inter- 
mediate shapes with bumps. We characterize these new 
equilibrium shapes in terms of order parameters (rota- 
tional invariants), which were introduced to describe virus 
capsids [2] and bond-orientational order in liquids and 
glasses |15j . Moreover, by comparing rotational invariants, 
we draw an analogy between the buckling transition found 
in virus capsids |14j and the transition from spherical 
shapes towards shapes with dimples and ridges explored 
in this paper. 



2 Phenomenological models 

The spontaneous-curvature model introduced by Helfrich 
in 1973 [8 accounts for a possible asymmetry of mem- 
branes, such as a difference in the number of molecules in 
each layer of bilayer vesicles. He suggested the following 
free energy of fluid membranes 



J r H = Jj dS {2k{H - H ) 2 + kK}, 



(3) 



where Ho is called the spontaneous curvature, k is known 
as bending rigidity and R is the Gaussian rigidity that 
affects only transitions implying a change of topology. 
The Helfrich model was successful in describing differ- 
ent phenomena, like the budding transition, discocyte- 
stomatocyte transition, et cetera |9I10| . However, for some 
systems it turns out to be too simple and thus insufficient 
to explain new experimental data. In fact the Helfrich 
approach accounts only for the basic degrees of freedom, 
such as bending, neglecting the possibility of a tilt of the 
molecules within a layer [16] , stretching/compression of 
layers |17) and interaction with the environment . The 
latter situation may result in a free energy with sponta- 
neous bending (negative coefficient in front of H 2 ) given 
by Eq. @. In general, one can imagine other mechanisms 
favouring spontaneous bending, like geometric constraints 
of packing, complex van der Waals or electrostatic inter- 
molccular interactions. 

In this paper we consider only symmetric fluid mem- 
branes, where the Helfrich free energy J-~h with Ho = 
coincides with the Willmore functional W = jj dSH 2 . 
It is worth noting that J-"h does not depend on the size 
of the structure. The presence of higher order terms in 
J-4 (Eq. ([2])), on the contrary, result in a characteristic 
length scale, L oc ^Jci/A. To guarantee that T4 has a 



minimum for real values of H and K we require the form 
${H 2 ,K) = ciH A + c 2 H 2 K + c 3 K 2 to be positive definite 
for H 2 ^ K . Thus a minimum exists if 



ci > 0, c 3 > 0, 4ciC3 



Co > 0. 



(4) 



To relate the coefficients of Eq. @ to the bending rigid- 
ity k, entering the Helfrich model Eq. Q, we calculate 
d 2 Fi/dH 2 at H = Hq which provides a minimum of the 
functional. The result is 



4:Aqc 3 
4ciC3 — a 



> 0. 



(5) 



In the case c 2 = 0, Vci, C3 > this expression simplifies to 
K = A. In principle the parameters Ci can be determined 
by comparing equilibrium shapes with experimental val- 
ues for the deformations due to high magnetic fields as 
it has been done for sexithiophene vesicles |12j . However, 
the coefficients Ci are not separately accessible. Here, we 
will consider them as formal parameters, satisfying Eq. Q 
and study the possible equilibrium shapes and their trans- 
formations, without referring to specific systems. 

To find equilibrium shapes, one usually needs to solve 
the Euler-Lagrange equation. For the Willmore functional 
W only six analytic solutions are known: planes, cylinders, 
spheres, tori, cones, dupin cyclides [IS]- Since we do not 
consider here topological transformations, the sphere gives 
the absolute minimum W = 4-7T for x = 2. The equilib- 
rium shape equation for the functional J-4 can be written 
following Ref. [TS], which results in 

2AH{H 2 -K)+A{c 2 -c l )H z K+2{c 3 -c 2 )HK 2 +Qc l H b + 

+ V 2 {AH + 2qH 3 + c 2 HK) + W 2 (c 2 H 2 + 2c 3 K) = 0. 

(6) 

Here V 2 stands for the Laplace-Beltrami operator, V 2 = 
(^/ \/\9\)9i(g^ \/\g\di) , gij is the metric tensor of the sur- 
face,^' = (g kl )7.\ \g\ = det \\ gij \l V 2 = (l/VR)Si(XL« 



|<?|<9j), Lij is a second fundamental form [B]. By sub- 
stituting H 2 = K = l/R 2 , where R is the radius of a 
sphere, into Eq. ([6]) , one finds a sphere as solution if and 
only if ci + c 2 + c 3 = 0, which contradicts the condi- 
tion (j?]). Therefore a sphere is not an equilibrium shape 
of the free energy ([2]). Finding analytical solutions for a 
highly nonlinear partial differential equation (Eq. ©) is 
not likely. Alternatively, one can study equilibrium shapes 
by means of numerical methods. Here we adapt the "Sur- 
face Evolver" software [13120] to this purpose. 



3 Computational details 

In the interactive program "Surface Evolver" , a surface is 
modeled by a set of triangles, which is finite for compact 
surfaces. Given a triangulation of the surface we evolve it 
towards the shape that minimizes the total free energy. For 
the evolution "Surface Evolver" uses the steepest descent 
method, which means that at each iteration all vertices 
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Fig. 1. An example of the evolution of icosahedron (a) towards a sphere (c) for the Willmore free energy. Every 
iteration, each triangle is subdivided into four similar smaller triangles. The blue color on (c) illustrates the vertices 
with a fivefold coordination. 



are moved along the gradients of the free energy. Then we 
can refine the surface by dividing each triangle into four, 
and repeat the procedure until the approximated surface 
becomes smooth. An example of the evolution, starting 
from the icosahedron, towards a sphere is shown in Fig. [TJ 
Since the Euler characteristic x does not depend on the 
triangulation of the surface, but only on the topology, it 
always holds \ — E + V = 2, which relates the total 
number of triangles (faces F), edges (E) and vertices (V). 

The first application of "Surface Evolver" was to study 
the shape minimizers of the Willmore functional W, start- 
ing from polyhedra with different \ [4] ■ The Helfrich spon- 
taneous curvature model (Eq. ©) with Hq ^ and vol- 
ume constraint was analyzed with this program in Ref. |21) , 
resulting in the prediction of corniculate, tubelike and 
other nonaxisymmetric shapes of vesicles. Among other 
examples are the study of the rheology of foams, simu- 
lation of microgravity, phenomena of capillarity and wet- 
ting [20] . 

Here we aim at studying the energy minimizing shapes 
of the free energy ([2]), starting, for reasons that we explain 
later, from an icosahedron (see Fig. la). Surface Evolver 
v.2.30 evaluates the energy terms, JfdS H 2 and JfdS K 2 . 
We have written two new subroutines to calculate the 
other two terms entering Eq. (J2J), namely fj dS H 4 and 
J J dS H 2 K [22]. For discrete surfaces at each vertex v, 
the mean curvature H u and the Gaussian curvature K u 
are defined as [23124) 

Bu= ] ^, K^^Hk-^M, (7) 

i 

where S v is a Voronoi area around v, i.e. the area of the 
cell consisting of all points closer to v than to any other 
vertex, VSV is the gradient of S v at v, and J^. is the 
sum of all facet angles at the vertex v of icosahedron. The 
definition of H v comes from the fact that the mean (extrin- 
sic) curvature measures the variation of the area element, 
displaced along the normal, divided by the correspond- 
ing change of volume. The definition of K v comes from 
the Gauss-Bonnet theorem for a Voronoi region. Then, 
we can approximate the integrals by assigning their en- 
ergy contributions to the vertices only. The integrals are 
calculated as a sum over all vertices v of the curvature 



times the area around a vertex, which gives: 




(9) 

Note that, for the calculation of the energy contributions, 
the area associated with a vertex is taken to be 1 /3 of S v , 
rather than S v . This approximation simplifies the calcu- 
lations for "Surface Evolver" and works best for triangles 
close to equilateral. The choice of icosahedron as a starting 
shape guarantees that, at every iteration, we are close to 
equilateral triangulation so that the approximation of the 
integrals as in Eqs. (JSJ and (0) holds. Moreover, among 
all regular polyhedra, the icosahedron has the ratio of the 
surface area over the enclosed volume closest to that of a 
sphere. It is in general convenient to start the minimiza- 
tion from a shape close to a sphere because the steepest 
descent method implemented in "Surface Evolver" finds 
only one stable local minimum, which is not necessarily a 
global one. Moreover, in nature many shapes, like viruses, 
are found to have an icosahedral symmetry and we will 
compare the predictions of our phenomenological model 
with the models developed for viral capsids, in particular 
the one discussed in [T4"] . 



4 Numerical results 

Assuming a constant area of a surface S = An (R = 1), 
we minimize the free energy given by Eq. ([2J with "Surface 
Evolver" . In Table Q] we illustrate the equilibrium shapes 
as a function of A for three typical cases: i) c\ = 1, 
C2 = C3 = 0, ii) ci = 1, C2 = 0, and C3 = 0.5, iii) 
ci = 1, C2 = — 1, C3 = 0.5. Irrespective of the partic- 
ular choice of the constants spherical shapes (shapes 
very close to a sphere [25]) are found for negative values 
of A (~A J J dS H 2 is positive). For A~4we see a tran- 
sition towards dimpled shapes (second column in Table [IJ 
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Table 1. Summary of the equilibrium shapes for different values of the parameters in T/± = JJdS{—AH 2 + ciH 4 + 
C2H 2 K + C3K 2 }. In the rows from left to right, shapes evolve from spherical towards non-spherical ones with increase 
of the negative contribution —A J J dS H 2 . The number of the column, given in brackets, corresponds to the label 
of the points in Fig. [3] and Fig. 2] All the shapes have a constant area and triangulation with a number of vertices 
V = 642. 

(point) 



c 2 = 
c 3 = 



c 2 = 
c 3 = 0.5 



C2 
C3 



-1 

0.5 




followed by a transition towards shapes with ridges (last 
column in Table [IJ connecting the vertices of icosahedron. 
In presence of a non-zero coupling term ci between the 
mean curvature H and the Gaussian curvature K, inter- 
mediate shapes with bumps (see bottom row in Table [T]) 
occur, favouring a positive K for C2 < 0. Note that, the 
bending rigidity n increases with A according to Eq. 
Thus the stiffness of shapes increases in the rows, for the 
first and the second rows n — A (02 = 0), for the third 
row k — 2A. 

A canonical way to characterize the shapes presented 
in Tablc[T]is to expand their radial density R(0, <f>) in terms 
of spherical harmonics Yi m (9, </>), as follows, 



R(0, 



(10) 



1=0 : 



-—I 



with the coefficients Qi m of the above expansion defined 
as 



Qi 



dej) / d9 sm6 R(6,(f>)Y l l l (6,<p). 
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In the case of triangulated surfaces we deal with a discrete 
radial density R(9, 4>), defined at the vertices of triangles. 
The coefficients Qi m were computed in Matlab using Gauss 
quadratures. Then, according to [15], we construct second 
and third order rotational invariants, as 



Qi 



-in 



1 \{2l + l) 



(12) 



and 



Wll ,7112 ,7713 [ rt S rt ) Qlm t Q hn 2 Qlm 3 
mi +7712+7713=0 V TOl TO2 



3/2 



(13) 



where 



are the Wigner 3j symbols. For a 



I I I 

mi m,2 TO3 

sphere the only non-zero coefficient is Qoo = V47r, giv- 
ing Qo = 47T. The shapes with icosahedral symmetry are 
distinguished by non-zero invariants Qi and Wi for I = 
0, 6, 10, 12, . . . [T3]. The vanishing of invariants for other 
values of I was also recovered in the present calculations. 
For the icosahedron, only for very high refinements, namely 
for triangulation with number of vertices V = 10242, we 
found for Wq and W\q (notice Wo = 1) the same values as 
in Ref. [15]. Therefore, in the following, all the integrals in 
Eq. (fTTj) are computed with V = 10242, rather than the 
one presented in Table [1] 

We plot first Qo/Att, characterizing the aspherity of 
the shape, as a function of the adimensional parameter 
AR 2 jc\ in Fig. [5] We define spherical shapes [3S] when 
Qo/Att ~ 1, implying that AR 2 fc\ ^ 2. Non-spherical 
shapes correspond to AR 2 jc\ > 2 for all three curves with 
different values of C2 and C3. The biggest deviation from 
a sphere, Qo/Att m 0.86, happens at A = —12, c 2 = and 
C3 = 0.5. The corresponding dimpled shape can be found 
in Table [1] Notice, that the volume of the equilibrium 
shapes decreases with increasing A in the same way as Qo ■ 
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Fig. 2. Aspherity of the shape, described by the lowest 
spherical harmonic Qq, normalized to the value Qq — 4ir 
of a perfect sphere, as a function of AR 2 jc\ with c\ = 1, 
R = 1 (see Eq. [5]) . The vertical line indicates the crossover 
between spherical and non-spherical shapes. The value of 
Qo/Att ~ 0.967 for the icosahedron is shown as horizontal 
line. 




Q 6 /Qo 



Fig. 3. The second order invariants defined in Eq. (TT2"j) . 
The points 1-6 correspond to the values of A — 
0,4,8,12,16,20 respectively, as in Table Q] For compar- 
ison, we show the data from [Tl], used to describe the 
buckling (points b-d) and faceting (points d-f ) transitions 
of viral capsids. 



Fig. 4. The third order invariants We and Wiq distinguish 
the shapes with bumps and ridges (We < 0) from the 
dimpled shapes (We > 0). The points 1-6 are labeled as 
for Fig. |3] and Table [TJ The spherical shapes correspond to 
point 1 at (0, 0). The icosahedron is marked by an asterix. 

less Foppl-von Karman (FvK) number 7 = YR 2 /k, where 
Y is the Young modulus. The points b-d describe the 
buckling transition of viral capsids, and the points d-f are 
associated with sharpening of the ridges at large 7 [14] . 
In our case, the transition towards dimpled shapes may 
be analogous to the 'buckling', whereas the appearance 
of ridges leads to the increase of Qw/Qo contrary to the 
model of viral capsids. 

To find out more connections between the non-spherical 
equilibrium shapes and to distinguish shapes with bumps 
and dimples, we consider the combination involving the 
third order invariants W\ , as shown in Fig. 0] We notice 
that shapes with dimples are characterized by WeQe ~ 
0.02, whereas shapes with bumps have the same magni- 
tude but the opposite sign of We (curve with rhombuses). 
The sign of We, which is the first non-zero third-order in- 
variant, discriminates the shapes with dimples from the 
shapes with bumps and ridges. The sharpening of ridges 
is characterized by an increase in the invariants Qi and Wi 
calculated for higher degree Z = 10. According to Figs. OS] 
the point 3, corresponding to A = 8, c\ = 1, c% = — 1, C3 = 
0.5 (see also Table]!]) is the closest one to the icosahedron. 



Figure [3] shows the second order invariants Q1/Q0 for 
I = 6, 10 plotted against each other. Starting from a spher- 
ical shape labeled by point (1) we follow the curves with 
the points separated by S A — 4, as in Tabled] The dimpled 
shapes correspond to the growth of Qe / Qo (point 2 for all 
curves) , whereas the appearance of ridges is characterized 
by the increase of Q10/Q0 along the curves. By adding 
the points from Fig. 8 in Ref. [14] , we compare our equi- 
librium shapes for fluid membranes, with shapes of crys- 
talline viral capsids studied within the continuum elastic 
theory. The governing parameter of that model, relating 
stretching and bending of elastic shells, is the dimension- 



5 Discussion and conclusions 

We have studied the equilibrium shapes of symmetric fluid 
membranes with a spherical topology. Assuming a fourth- 
order curvature model proposed in Ref. |11] we found a va- 
riety of shapes with dimples, bumps and ridges as well as 
quasi-spherical shapes. We noticed that similar shapes ap- 
pear in the theory of elastic icosahedral shells, when study- 
ing the buckling and ridge-sharpening transitions |14l26j . 
Both these transitions depend only on the FvK number, 
which is the ratio of the stretching and bending contri- 
butions to the free energy. In our case, the competition 
arises between the negative quadratic term in Eq. ^ 
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and higher order quartic terms. Our numerical analysis 
shows that the transition from spherical towards dimpled 
shapes depends only on the value of A or more likely on 
the dimensionless combination AR 2 jc\. The transitions 
between non-spherical shapes, such as dimples-ridges and 
dimples-bumps, are not determined only by the parameter 
AR 2 1 c\. Both buckling and ridge-sharpening transitions 
occur within one order of AR 2 j c\ whereas in the model of 
elastic shells the FvK number should change within four 
orders of magnitude [14126] . 

Our calculations were done under the constraint of 
constant surface, and we found a decrease of the volume 
with increasing A, similar to the change of volume upon 
buckling transition [26 . It might be interesting to study 
the ^-minimizing shapes under the constraint of constant 
volume. With the latter constraint, non-spherical shapes 
appear also with the Willmore functional but those shapes 
are essentially different from the ones we find within our 
model, namely they present large deviations from spheres 
but no bumps, ridges or dimples. 



20. Http://www.susqu.edu/brakke/evolver/evolver.html 
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